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, Abstract 

^} ■ We apply the iterative numerical algorithm ALPHA, which automatically gen- 

erates the N-point Green functions (as recently suggested Q) to the computation 
of the production rate of four fermions and four fermions plus photon in electron- 
positron annihilation. The discussion of the physical results is preceded by an 
introduction on the algorithm. 
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1 Introduction 

The center of mass energy available at the LEP I collider will be well above the threshold 
energy for W pairs productions and one of the main goals of the collider is to allow the 
study of the W physics (as well as of any new physical effect which might manifest at the 
200 GeV scale). The expected integrated luminosity is about 500 pb^ 1 per experiment and 
this will allow a relatively good measurement of the W pair cross section and therefore 
of the weak boson mass. The W decays into two fermions and thus the signal for such 
a process is the presence of four fermions with the appropriate charge (and flavor). The 
narrow width approximation namely the assumption that the rate is given by the rate for 
W pair production times the branching ratio of the W decay receives important corrections 
already at the tree level for two main reasons: i) the production is close to the threshold 
and the presence of a lineshape for the W propagator significantly affects the rate, ii) other 
graphs contribute to the same final state in addition to those accounting for virtual W 
exchange. Both correction are large enough to be 'detectable' with the expected statistics 
and therefore a full computation is required. 

*INFN Fellow 
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We also present the calculation of the production of four fermion plus a detectable 
photon arising from anyone of the charged virtual/real particle involved in the process. 

To perform the computation of the cross sections we have used the code ALPHA |l[ 
for the scattering matrix elements and the package VEGAS |2[ for the integration over 
the phase space. 



1.1 The Method for the Computation of the N-points Green 
Functions (ALPHA) 

In a recent paper, |l[ , we have suggested an iterative algorithm to compute automatically 
the scattering matrix elements of any given effective lagrangian, T. The S matrix is re- 
covered from the effective lagrangian after a functional legendre transform. By exploiting 
the relation between T and the connected Green function generator, Z, we have obtained 
an iterative formula which does not require the explicit use of the Feynman graphs and 
it is suitable to be implemented in a numerical routine. The problem of computing a 
N-point scattering matrix element can be recast as the problem of finding the minimum 
of a function Z with respect to a finite set of variables. Once the subsequent stationary 
conditions for Z are written down they can be solved iteratively (with respect to an ex- 
pansion parameter). After a proper number of iteration steps one obtains the scattering 
amplitude. 

The objects of interest to us are the connected Green Functions of the theory Zj^) =< 
4>(j>\) ■ ■ -<j)(pn) >• To compute one of them, it is convenient to define a function Z from 
the functional Z . By means of a particularly simple choice of the source term (following 
the notations in [l|]) 

J(x) = a ie -^ x (1) 



we define 



/n 
dxT((j)(x)) + ^2 a j ( t ) {Pj 



sz 

(3) 

where the stationary condition (|2|) implies that (f) has to be considered as a function of 
the a,j such to satisfy (13): in other words if the field <f>(x) is written in terms of its fourier 
components 

<j>{x) = b ie - lP ' x (4) 

then the bi are function of the since they satisfy the equations dZ/dbi = 0. The index 
% runs over a discrete set of values. In fact, at the stationary point (|2]), only the momenta 
Pi which are a linear combination of the pi with integer coefficients have a non zero hi 
(see after). Moreover the Green Function Z^l, pn is immediately found as 
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and only a finite number of bi is relevant in this limit. Indeed the bi can be expanded in 
a powers series in the a*, 

oo 

h= E i>r"- L "'r. (6) 

h=0,-,ln=0 

and all the terms containing a™ (n > 1) do not contribute to the limit (||). In practice for 
a lagrangian of the type 

£ = -^^fl^ + -^^^0^ (7) 
2 6 

the stationary condition (0) will take the form 

< = fl£6& + ^V^Z- (8) 

where we have denoted with greek index all possible particles degrees of freedom: lorentz, 
color, . . . This equation can be solved perturbatively with respect to the a, or, equivalently, 
to the interaction terms O: if we call b^ m the solution up to order m, the following recursive 
relation holdsQ 



bf fi = a" j = l,n 

6J = j > n n = number of external particle 
b a --Y\ af3 n^ 5 V b 1 b s 

u j,m — 9 t,k,l u k,r u l,s 

~ r+s=m— 1 

(9) 

and the amplitude is given byf] 

z s+r=n-2 s+r+t=n-3 

(10) 

Each has a precise physical meaning: it represents the sub-amplitude with m vertex 
interactions and with the following external legs: one external leg with off-shell momentum 
Pi plus any subset of external on-shell momenta chosen among the external momenta p^. 
The bj tm is zero if there is no subdiagram with m vertex and with the required external 
legs. 

Thus the above algorithm is a recursive relation between the subdiagrams of the full 
process: all the off-shell subamplitudes with m vertex are given in terms of all the sub- 
amplitude with k < m vertices. 

the first step we have simplified a factor f[°^ to obtain the truncated Green function. 
2 One actually needs to know a number of b" lower than the one which can be inferred from eqs. (^|) 
and (^|) since the stationary condition can be used to further reduce the number of relevant variables. 
We refer the interested reader to pi. 
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It is manifest that this recursive formula does not require any symbolic manipulation 
and does not require the knowledge of the Feynman diagrams. Once the numerical values 
of the external momenta pj are specified the b^ r and O^jfj in (§) and ((TOj) are just arrays 
of complex numbers and therefore this formula can be used to implement a numerical 
code. 

1.1.1 The Code ALPHA 

When the initial and final states of the process are specified (type, momenta and spin 
of the external particles) the program ALPHA prepares an array 6" for all the possible 
degrees of freedom. 

The coefficients and LT^ are returned by some subroutines as a function of the 

finite set of possible momenta P m . 

The ALPHA code includes all the electroweak interactions and all the flavor content of 
the Standard Model (SM) and it can perform the computation of all electroweak matrix 
elements in the SM regardless of the initial or final state type. In addition, due to its simple 
logic, it allows to modify the lagrangian with no excessive effort (by adding the proper 
subroutines to compute the new interactions and/or adding the relevant variables for 
the new particles). 

The algorithm defined from the recursive relations ([U]) is fast and the computing time 
of the code seems to increase roughly as a n (with a ~ 4), n being the number of external 
particles, while the number of feynman graphs (and therefore the time for the evaluation 
of their sum) usually grows faster than n\. 

The origin of this difference might be understood by comparing the number of opera- 
tions required in the evaluation of the left-hand side and right-hand side of the equation 
a{pi + ■ ■ • + b n ) = ab\ + • • • ab n ; to compute the sums before the multiplications is more 
economical than the inverse. On the contrary one uses to compute all the Feynman graphs 
and only at the very end sums all their contributions. However, especially for complicated 
processes, the same sub-amplitude often occurs in several distinct Feynman graphs. In 
the above algorithm this factorization is explicitly exploited since the b^ m do represent 
such sub-amplitudes. 

Finally, since the code is entirely numerical, the output can be immediately used for 
the integration procedure. We will comment in the following sections about the reliability 
of the code and its performances. 

1.2 The integration procedure 

In order to construct a kinematically allowed momenta configuration the phase space is 
factorized as a multiple decay process using the formula || 

d 3 "" 4 $(P; q u q 2 , q 3 , q n ) = d 2 §{Q = gi + q 2 ; q x , q 2 )d 3n - 7 <5>(P; Q, q 3 , q n ){2ixfdQ 2 (11) 

whose main advantage is that one just need a single subroutine for a two body decay and 
then by repeated calls to this subroutine it is straightforward to obtain an N body phase 
space. 
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To compute the cross sections of interest one must overcome some problems Q. When 
the number of phase space variables is large the only possible approach to the phase space 
integration is via MONTECARLO integration. The scattering amplitudes can be singular, 
or more generally strongly peaked, in particular regions of the phase space. This can be 
due, for example, to the presence of infrared (quasi-) singularities, or to the propagation of 
heavy unstable particles. When this is the case some of the internal propagators become 
huge and, as a result, the bulk of the cross section comes from a very thin region of the 
phase space. To deal with this problem it is mandatory to have a large sample of points 
in the singular regions. To this purpose we have used the package VEGAS [0 whose 
underlying idea is to redefine the integration variables in order to minimize the variance. 
So, for example, for the process e + e~ — > h + u^t^u t we have chosen as phase space variables 
the invariant masses \L\ and \ii of the fermion pairs and t~v T) the angle 9 between 
the beam and the combined t~v t momentum plus other five variables, along the lines of 
(PI), which don't exhibit peaking behavior. VEGAS then rescales properly /i 2 and 9 
to account for the peaking behavior of the cross section, namely it greatly increases the 
sampling for \i\ and /i 2 close to the mass of the W boson and it increases the sampling 
close to 9 = 0. The situation is more delicate when the number of singularities is bigger, 
mainly when there are many almost on-shell virtual particles which can lead to the same 
final state. This is the case, for example, for the process e + e~ — ► /i + z/^/i~z/ M which can 
proceed both via W + W~ and via Z Z or Z 7* virtual production and decay. VEGAS is 
efficient when the '(quasi-) singular' piece of the amplitudes can be written as the product 
of functions fj{xj) (where dfj/dxk ~ Sjk) of the integration variables Xj, thus a simple 
choice of phase space variables along the lines of (|TT|) is not satisfactory. To circumvent 
the problem we have split the integration domain into two regions and we have performed 
the integral with two different choices of phase space variables namely: 

a) in the region where virtual W exchange is dominant we have used: the invariant 
masses \i\ and \x% of the fermion pairs [i^v^ and /i - ^ the angle 9 between the beam and 
the combined [x'v^ momentum plus other five variables, 

b) in the region where virtual Z or 7 exchange is significant we have used: the invariant 
masses ii\ and /i 2 of the fermion pairs yU + /i~ and v^v^ the angle 9 between the beam and 
the combined fi~n + momentum plus other five variables. 

To precisely define the two distinct regions we evaluate the inverse propagators of the 
virtual W + , W~, Z and 7 and we use the variables in a when the product of the W 
propagators is smaller than the one of the Z or of Z and 7, otherwise we use the variables 
defined in b. 

An analogous problem occurs for the process e + e _ — > e + u e fi~u fl . In this case in 
addition to the usual W pairs exchange an additional peak in the cross section occurs when 
the final e + is emitted in the forward direction with respect to initial one, thus allowing 
for a soft virtual photon to be emitted. Again we split the integration domain into two 
regions defined by the relative size of the virtual W + and of the virtual 7 (invariant e + , 
e + four momentum) propagator and we perform the integral using two different phase 
spaces. 

It is straightforward (only a bit tedious) to extend the procedure to final states con- 
taining e + e~ pairs with the obvious drawbacks that the integration domain needs to be 
split into a larger number of regions. 
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The processes with the emission of an additional final photon have additional singu- 
larities. It is well known that they are dominated by the emission of soft collinear photon 
and this remain true even with the relatively hard cuts we impose on the photon four 
momenta. One needs therefore a phase space which allows to single out the soft and 
collinear behavior of the cross section. Let us consider the process e + e~ — > T + z/ 7 -/i~z/ / /y. 
The phase space is divided into three regions according to which among the invariant 
masses of each charged fermions with the photon is the smaller. If it happens to be the 
e + 7 or e~7 invariant mass than we choose as integration variables the momentum of the 
photon its angle with respect to the beam and the invariant W + , W~ masses and the 
angle of the W with the beam. If the photon is more collinear with the r + we choose as 
variables the invariant masses [i~v^ r + i/ T 7 and t + 7, then the momentum of the photon 
and the angle between the W~ and the beam. Analogously if the photon is emitted closer 
to the n~ (these two last region are absolutely symmetric and the same VEGAS' sampling 
can be used). 

The procedure in other cases is completely similar. 

2 e + e~ Into Four Fermions (plus 7) 

Before computing the four fermion processes (plus 7) we have to fix the interactions in 
the lagrangian and their couplings. We consider the Standard Model lagrangian; the 
interactions involved in our processes include all the gauge bosons three-vertices, some 
four-vertices (for 4/ + 7) and all these interactions are unambiguously defined once we 
require them to satisfy the 577(2) x U(l) symmetry conditions and we specify the couplings 
g>2 = 0.651698, gi = 0.357213 and the mixing angle angle 9w = arctan(gi/ g-i) between 
the two neutral gauge boson. The couplings among the gauge bosons and the fermions 
are chosen as follows: the photon fermion interaction is settled to a(2My/) = e 2 /(47r) = 
1/128.07, the VT-fermion interaction is given by the SU(2) coupling g%. Finally the Z- 
fermion vector gy and axial- vector g& coupling are derived from the SM relation gv/ gA = 
1 — 4\Q\sin 2 9w (Q is the charge of the fermion in unit of electron charge) and gA = 92/^- 
Gluons are not taken into account. The gauge boson masses (in GeV) are Mz = 91.1888 
and My/ = 80.23; we added an imaginary part to the masses in the propagators in order 
to take into account of the constant (or running) || width of the physical weak boson 
propagator (Tw — 2.4974, Y z = 2.03367). All fermion masses are taken from the Review 
of Particle Properties p3] . 

2.1 Results and Comments 

Whenever possible the four fermion processes are dominated by the production of a W 
pair with the subsequent decay. The naive narrow width approximation (i.e. to replace 
the cross section by the cross section of the process e + e~ — > W + W~ times the branching 
ratio of the W into the given channel) receives important corrections already at the tree 
level due to the finite W width and to the presence of non resonant diagrams leading to 
the same final state as those mediated by the W decay. The effect of the finite width 
is typically of several per cent and the effect of non resonant diagrams is particularly 
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important for processes where the number of non resonant diagrams is high and other 
channels contribute (Z and soft virtual photon) like for example e~z/ e /i + z/ M or yTv^^v^ 
(and similar semileptonic or hadronic channels). Moreover the non resonant diagrams 
become more important close to the threshold for W pair production namely into one of 
the energy regions which will be used for an accurate determination of the W boson mass. 
Because of the high accuracy aimed to at the LEP II experiments all these effects need 
to be accounted for. In table [I] we list the tree level cross sections for several four fermion 
final states at a cm. energy of 190 GeV. Within the Montecarlo error the results are in 
agreement with those obtained by means of other methods [6] thus proving the reliability 
of our code and its usefulness to build up Montecarlo Event Generators for complicated 
processes. It is important to stress here that the whole procedure for the computation 
of the matrix element is completely automatic; once the algorithm is implemented (like 
in ALPHA for example) the user has only to specify the number and type of external 
particles in an input file. As output the code will directly give the scattering matrix 
elements for any given spin and momenta configuration. In table |3| we report the CPU 
times required for the event generation for ALPHA to assess the performance of our code 
(although one should keep in mind that it is just the first implementation of the algorithm 
and large room for optimization is left). 

The addition of a photon in the final state introduces a factor oqed m the total cross 
section, but the existence of infrared and collinear singularities, allows large logarithms 
to appear and to give a sizeable cross section. 

The most important diagrams to contribute are those accounting for the emission of 
the photon from one of the external charged fermions. In such a case the amplitude gets 
a factor 

■jzi — * U W ^ 

f + If —m 

coming from the internal propagator contiguous to the external photon. When (p + k) 2 — 
m 2 goes to zero, the amplitude increases. This introduces into the cross section a a peak 
behavior with respect the angle 9 (between the photon and the charged fermion) and the 
energy E n of the outcoming photon (neglecting suppressed terms of the type m/p ■ k) 

dE^d cos # 7 

* ~ 25(1 - cos ft,)" (13) 

We have introduced a cut-off to keep under control these singularities and to take into 
account the experimental detector capabilities. If we call d{ the angle between the photon 
and the external charged % - th fermion, then for each final state configuration, we define 
6ty as the minimum angle among the 

In figures (1-3) we plot the differential cross section for different processes as func- 
tion of £7 7 and # 7 after having integrated over all the other phase space variables. The 
corresponding cross sections are reported in table |2|. We do report only a subset of the 
cross-sections for e + e~ into four fermion plus 7. The given processes, however, exhibit 
the most general peaking structure for this class of processes (with the exception of those 
containing an e + e~ pair into the final state) and therefore the extension of the calculation 
to any other process of this type is immediate. 
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During the LEP 2 phase these processes will be observable and the the hard photon 
emission could be measured. We have checked that the overwhelming contribution to 
these cross sections is due to the emission of soft and/or collinear photons and that the 
four point gauge boson vertices are negligible. If the born cross section is multiplied by 
the eikonal factor (and the proper rescaling of couplings) one recovers with a relatively 
good approximation (from few to ten percent in the relevant regions) the shown plots. 

In table |3] we report the CPU performances of the code for a few four fermion plus 
7 final states. They are typically a factor three larger than the corresponding ones for 
the processes without photon emission. In view of the large increase in the number and 
in the complexity of the Feynman graphs we regard this as remarkable property of our 
algorithm. 
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process 


cross section (pb) 


process 


cross section (pb) 


e+e~ — > u, + v,,t~Vt 


0.2185(1) 


e+e" — > ii^VnUd 


0.6663(4) 


e+e" — ► udsc 


2.015(1) 


e+e" — > e _ i/ P iA,u+ 


0.2271(1) 


e+e" — > e~v e ud 


0.6929(5) 


e+e - — ► u + u VuVii 


0.2286(2) 


e+e" — > uudd 


2.063(1) 


e+e" — ► e+e"z/ e z/ e 


0.2573(2) 


e+e — — > /i+ ji~ v T v T 


0.01010(1) 


e+e - — >■ fi + fi~T^T~ 


0.00925(1) 
0.01100(1) 


e + e~ -> U^U^rUr 


0.008245(4) 


e+e~ — > n + n~uu 


0.02448(3) 


e+e~ — > fi + fi~dd 


0.02374(2) 


e+e" — > v^v^uu 


0.02104(2) 


e+e" — > v^v^dd 


0.01988(2) 


e+e" — > wucc 


0.05221(5) 


e+e" — > muss 


0.04958(5) 


e+e~ — ► oMss 


0.04704(5) 


e+e" — > u e u e fi + fi~ 


0.01778(2) 




0.008335(4) 


e+e~ — > v e v e uu 


0.02389(2) 


e+e~ — * u e u e dd 


0.02066(2) 


e+e" — > fj, + fj,~ fi + /i~ 


0.005456(5) 


e+e" -»■ v^v^v^v^ 


0.004065(5) 


e+e" — > wiZim 


0.02565(3) 


e+e" — ► dddd 


0.02349(2) 


e+e" - ► v e v e v e v e 


0.004091(2) 





Table 1: Cross sections in pb for various 4 fermions final states in e+e~ collisions at a 
center of mass energy of 190 GeV. The following cuts are applied: The invariant mass of 
each qq, qq and qq pair is greater than 5 GeV, the energy of each final quark or antiquark 
is greater than 3 GeV, the energy of each final charged lepton is greater than 1 GeV, the 
angle of each charged leptons with the beam direction is greater than 5°. The standard 
model couplings are defined in the test and all fermion are massive with the masses as in 
[3]. For the final state /i+/x"r+T" the cross section is given for both massive and massless 
fermion (lower row). In all other cases the difference is small. Gluon exchange is not 
accounted for. 



final state 






22.652(46) 


e + v e u~d~ r y 


58.00(17) 


u + d + c~s--i 


158.38(36) 


u + d + u~d~ r y 


163.43(47) 



Table 2: four-fermion + 7 final states in e+e" collisions. Inputs are the same as in table 
(1). The cuts are the same as in table (1), the photon energy is greater than 1 GeV, 
and the angle of the photon with the beam is greater than 11.5° Gluon exchange is not 
accounted for. 
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final state 


CPU, sec/100000 events 




270 




760 


e + v e ii~d~ 


445 


e + v e u~ d~ 7 


1388 


u + d + c~ s~ 


327 


u + d + c~s~'y 


968 


u + u~d~d + 


605 


u + u~d~d + j 


1700 


e + e~e + e~ 


1120 


e + e~e + e~7 


3399 


e + e _ e + e"77 


10729 



Table 3: Required CPU time for the computation (single precision) of 100000 matrix 
elements for various processes. All the Higgs couplings are set to zero. The times are given 
in seconds and the calculations have been performed with a DIGITAL machine ALPHA 
3000/600 with 64M of memory. 



11 



Figure 1: differential cross section (fb/BIN) of e + v e u~d~'-f as a function of 0.42 < 
cos(9j) < 0.98 and 1 GeV < Ej < 19.5 GeV ( see the text). All inputs parameters 
as in Table 1 
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Figure 2: differential cross section (fb/BIN) of u + d + c~s~ r ) as a function of 0.42 < 
cos(9 y ) < 0.98 and 1 GeV < E 1 < 19.5 GeV ( see the text). All inputs parameters 
as in Table 1 
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Figure 3: differential cross section (fb/BIN) of u + d + u~d~ r y as a function of 0.42 < 
cos(9 y ) < 0.98 and 1 GeV < E 1 < 19.5 GeV ( see the text). All inputs parameters 
as in Table 1 
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